#include<iostream>
#include<time.h>
#include<stdlib.h>
#include "matrixmul.h"
using namespace std;
void cin_matrix(double ** A, int row, int col){
	time_t t;
	srand((unsigned) time(&t));
    for(int i = 0; i < row; i ++){
        for(int j = 0; j < col; j ++){
            A[i][j] = rand() % 10 + rand() / double(RAND_MAX);
        }
    }
}
void general_mul(double** A, double ** B, double ** C, int m, int n, int k){
    for(int i = 0; i < m; i ++){
        for(int j = 0; j < k; j ++){
            C[i][j] = 0;
            for(int r = 0; r < n; r++){
                C[i][j] += (A[i][r] * B[r][j]);
            }
        }
    }
}
void based_on_software_mul(double** A, double ** B, double ** C, int m, int n, int k){
    for(int i = 0; i < m; i += 4){
        for(int j = 0; j < k; j += 4){
            C[i + 0][j + 0] = 0;C[i + 0][j + 1] = 0;C[i + 0][j + 2] = 0;C[i + 0][j + 3] = 0;
			C[i + 1][j + 0] = 0;C[i + 1][j + 1] = 0;C[i + 1][j + 2] = 0;C[i + 1][j + 3] = 0;
			C[i + 2][j + 0] = 0;C[i + 2][j + 1] = 0;C[i + 2][j + 2] = 0;C[i + 2][j + 3] = 0;
        	C[i + 3][j + 0] = 0;C[i + 3][j + 1] = 0;C[i + 3][j + 2] = 0;C[i + 3][j + 3] = 0;
			for(int r = 0; r < n; r++){
				register double B_i0_r_reg = B[r][j+0];
				register double B_i1_r_reg = B[r][j+1];
				register double B_i2_r_reg = B[r][j+2];
				register double B_i3_r_reg = B[r][j+3];
				
				register double A_i0_r_reg = A[i][r];
                C[i + 0][j+0] += (A_i0_r_reg * B_i0_r_reg);
				C[i + 0][j+1] += (A_i0_r_reg * B_i1_r_reg);
				C[i + 0][j+2] += (A_i0_r_reg * B_i2_r_reg);
				C[i + 0][j+3] += (A_i0_r_reg * B_i3_r_reg);
				register double A_i1_r_reg = A[i + 1][r];
				C[i + 1][j+0] += (A_i1_r_reg * B_i0_r_reg);
				C[i + 1][j+1] += (A_i1_r_reg * B_i1_r_reg);
				C[i + 1][j+2] += (A_i1_r_reg * B_i2_r_reg);
				C[i + 1][j+3] += (A_i1_r_reg * B_i3_r_reg);
				register double A_i2_r_reg = A[i + 2][r];
				C[i + 2][j+0] += (A_i2_r_reg * B_i0_r_reg);
				C[i + 2][j+1] += (A_i2_r_reg * B_i1_r_reg);
				C[i + 2][j+2] += (A_i2_r_reg * B_i2_r_reg);
				C[i + 2][j+3] += (A_i2_r_reg * B_i3_r_reg);
				register double A_i3_r_reg = A[i + 3][r];
				C[i + 3][j+0] += (A_i3_r_reg * B_i0_r_reg);
				C[i + 3][j+1] += (A_i3_r_reg * B_i1_r_reg);
				C[i + 3][j+2] += (A_i3_r_reg * B_i2_r_reg);
				C[i + 3][j+3] += (A_i3_r_reg * B_i3_r_reg);
            }
        }
    }
}
void cout_matrix(double ** C, int row, int col){
    for(int i = 0; i < row; i ++){
        for(int j = 0;j < col ; j ++){
            cout<<C[i][j]<<" ";
        }
        cout<<endl;
    }
}
void Matrix_Sum(int N, double** MatrixA, double** MatrixB, double** Sum_Matrix){
	for (int i = 0; i < N; i++){
		for (int j = 0; j < N; j++){
			Sum_Matrix[i][j] = MatrixA[i][j] + MatrixB[i][j];
		}
	}
}
void Matrix_Sub(int N, double** MatrixA, double** MatrixB, double** Sub_Matrix){
	for (int i = 0; i < N; i++){
		for (int j = 0; j < N; j++){
			Sub_Matrix[i][j] = MatrixA[i][j] -MatrixB[i][j];
		}
	}
}
void Matrix_Mul(int N, double** MatrixA, double** MatrixB, double** Mul_Matrix){
	for (int i = 0; i < N; i++){
		for (int j = 0; j < N; j++){
			Mul_Matrix[i][j] = 0;
			for (int k = 0; k < N; k++){
				Mul_Matrix[i][j] = Mul_Matrix[i][j] + MatrixA[i][k] * MatrixB[k][j];
			}
		}
	}
}
void Strassen(int N, double** MatrixA, double** MatrixB, double** MatrixC){
	double** MatrixA11;
	double** MatrixA12;
	double** MatrixA21;
	double** MatrixA22;

	double** MatrixB11;
	double** MatrixB12;
	double** MatrixB21;
	double** MatrixB22;

	double** MatrixC11;
	double** MatrixC12;
	double** MatrixC21;
	double** MatrixC22;
	//初始化
	MatrixA11 = new double*[N/2];
	MatrixA12 = new double*[N/2];
	MatrixA21 = new double*[N/2];
	MatrixA22 = new double*[N/2];

	MatrixB11 = new double*[N/2];
	MatrixB12 = new double*[N/2];
	MatrixB21 = new double*[N/2];
	MatrixB22 = new double*[N/2];

	MatrixC11 = new double*[N/2];
	MatrixC12 = new double*[N/2];
	MatrixC21 = new double*[N/2];
	MatrixC22 = new double*[N/2];
	for (int i = 0; i < N/2; i++)
	{
		MatrixA11[i] = new double[N/2];
		MatrixA12[i] = new double[N / 2];
		MatrixA21[i] = new double[N / 2];
		MatrixA22[i] = new double[N / 2];

		MatrixB11[i] = new double[N / 2];
		MatrixB12[i] = new double[N / 2];
		MatrixB21[i] = new double[N / 2];
		MatrixB22[i] = new double[N / 2];
		
		MatrixC11[i] = new double[N / 2];
		MatrixC12[i] = new double[N / 2];
		MatrixC21[i] = new double[N / 2];
		MatrixC22[i] = new double[N / 2];
	}
	 //将大矩阵分割为4个小矩阵
	for (int i = 0; i < N / 2; i++){
		for (int j = 0; j < N / 2; j++){
			MatrixA11[i][j] = MatrixA[i][j];
			MatrixA12[i][j] = MatrixA[i][j + N / 2];
			MatrixA21[i][j] = MatrixA[i + N / 2][j];
			MatrixA22[i][j] = MatrixA[i + N / 2][j + N / 2];

			MatrixB11[i][j] = MatrixB[i][j];
			MatrixB12[i][j] = MatrixB[i][j + N / 2];
			MatrixB21[i][j] = MatrixB[i + N / 2][j];
			MatrixB22[i][j] = MatrixB[i + N / 2][j + N / 2];
		}
	}
   //辅助矩阵S
	double** MatrixS1=new double*[N/2];
	double** MatrixS2 = new double*[N/2];
	double** MatrixS3 = new double*[N/2];
	double** MatrixS4 = new double*[N / 2];
	double** MatrixS5 = new double*[N / 2];
	double** MatrixS6 = new double*[N / 2];
	double** MatrixS7 = new double*[N / 2];
	double** MatrixS8 = new double*[N / 2];
	double** MatrixS9 = new double*[N / 2];
	double** MatrixS10 = new double*[N / 2];
	
	for (int i = 0; i < N / 2; i++) 
	{
		MatrixS1[i] = new double[N / 2];
		MatrixS2[i] = new double[N / 2];
		MatrixS3[i] = new double[N / 2];
		MatrixS4[i] = new double[N / 2];
		MatrixS5[i] = new double[N / 2];
		MatrixS6[i] = new double[N / 2];
		MatrixS7[i] = new double[N / 2];
		MatrixS8[i] = new double[N / 2];
		MatrixS9[i] = new double[N / 2];
		MatrixS10[i] = new double[N / 2];
	}

	Matrix_Sub(N/2, MatrixB12, MatrixB22, MatrixS1);		//S1 = B12 - B22
	Matrix_Sum(N / 2, MatrixA11, MatrixA12, MatrixS2);		//S2 = A11 + A12
	Matrix_Sum(N / 2, MatrixA21, MatrixA22, MatrixS3);		//S3 = A21 + A22
	Matrix_Sub(N / 2, MatrixB21, MatrixB11, MatrixS4);		//S4 = B21 - B11
	Matrix_Sum(N / 2, MatrixA11, MatrixA22, MatrixS5);		//S5 = A11 + A22
	Matrix_Sum(N / 2, MatrixB11, MatrixB22, MatrixS6);		//S6 = B11 + B22
	Matrix_Sub(N / 2, MatrixA12, MatrixA22, MatrixS7);		//S7 = A12 - A22
	Matrix_Sum(N / 2, MatrixB21, MatrixB22, MatrixS8);		//S8 = B21 + B22
	Matrix_Sub(N / 2, MatrixA11, MatrixA21, MatrixS9);		//S9 = A11 - A21
	Matrix_Sum(N / 2, MatrixB11, MatrixB12, MatrixS10);		//S10 = B11 + B12

	//辅助矩阵M
	double** MatrixP1 = new double*[N / 2];
	double** MatrixP2 = new double*[N / 2];
	double** MatrixP3 = new double*[N / 2];
	double** MatrixP4 = new double*[N / 2];
	double** MatrixP5 = new double*[N / 2];
	double** MatrixP6 = new double*[N / 2];
	double** MatrixP7 = new double*[N / 2];

	for (int i = 0; i < N / 2; i++)
	{
		MatrixP1[i] = new double[N / 2];
		MatrixP2[i] = new double[N / 2];
		MatrixP3[i] = new double[N / 2];
		MatrixP4[i] = new double[N / 2];
		MatrixP5[i] = new double[N / 2];
		MatrixP6[i] = new double[N / 2];
		MatrixP7[i] = new double[N / 2];
	}
	Matrix_Mul(N / 2, MatrixA11, MatrixS1, MatrixP1);		//M1 = A11 • S1
	Matrix_Mul(N / 2, MatrixS2, MatrixB22, MatrixP2);		//M2 = S2 • B22
	Matrix_Mul(N / 2, MatrixS3, MatrixB11, MatrixP3);		//M3 = S3 • B11
	Matrix_Mul(N / 2, MatrixA22, MatrixS4, MatrixP4);		//M4 = A22 • S4
	Matrix_Mul(N / 2, MatrixS5, MatrixS6, MatrixP5);		//M5 = S5 • S6
	Matrix_Mul(N / 2, MatrixS7, MatrixS8, MatrixP6);		//M6 = S7 • S8
	Matrix_Mul(N / 2, MatrixS9, MatrixS10, MatrixP7);		//M7 = S9 • S10

	Matrix_Sum(N / 2, MatrixP5, MatrixP4, MatrixC11); 
	Matrix_Sub(N / 2, MatrixC11, MatrixP2, MatrixC11);
	Matrix_Sum(N / 2, MatrixC11, MatrixP6, MatrixC11);
	Matrix_Sum(N / 2, MatrixP1, MatrixP2, MatrixC12);
	Matrix_Sum(N / 2, MatrixP3, MatrixP4, MatrixC21);	
	Matrix_Sum(N / 2, MatrixP5, MatrixP1, MatrixC22);	
	Matrix_Sub(N / 2, MatrixC22, MatrixP3, MatrixC22);
	Matrix_Sub(N / 2, MatrixC22, MatrixP7, MatrixC22);
	//合并为C矩阵
	for (int i = 0; i < N / 2; i++){
		for (int j = 0; j < N / 2; j++){
			MatrixC[i][j] = MatrixC11[i][j];
			MatrixC[i][j+N/2] = MatrixC12[i][j];
			MatrixC[i+N/2][j] = MatrixC21[i][j];
			MatrixC[i+N/2][j+N/2] = MatrixC22[i][j];
		}
	}
}

